function dT = chebyshev_derivative(N, x, n)
% 计算0到N-1阶第一类切比雪夫多项式在x处的n阶导数
% 输入:
%   N - 切比雪夫多项式的个数（计算T0到T_{N-1}）
%   x - 计算导数的点，应在[-1, 1]区间内
%   n - 导数的阶数（非负整数）
% 输出:
%   dT - N×1的列向量，包含各阶导数 T0^{(n)}(x), T1^{(n)}(x), ..., T_{N-1}^{(n)}(x)

% 初始化结果向量
%dT = zeros(N, 1);
dT = zeros(1,N);
% 检查n是否为非负整数
if n < 0
    error('导数的阶数n必须为非负整数。');
end
% 保证计算点在范围内
if abs(x) > 1+eps(max(1, abs(x))) * 10
    error('计算点x必须在[-1, 1]区间内。');
end


if n == 0
    % 当n=0时，直接返回切比雪夫多项式的值
    if N == 0
        return;
    end
    T = zeros(N, 1);
    T(1) = 1; % T0(x) = 1
    if N >= 2
        T(2) = x; % T1(x) = x
        for k = 3:N
            % 递推公式: Tk(x) = 2x*Tk-1(x) - Tk-2(x)
            T(k) = 2 * x * T(k-1) - T(k-2);
        end
    end
    dT = T;
else
    % 处理n>0的情况
    if N <= 0
        return;
    end
    
    % 计算最大的m = k - n，其中k的范围为0到N-1
    m_max = (N-1) - n;
    if m_max < 0
        % 所有k < n，导数均为0，直接返回
        return;
    end
    
    % 计算第二类切比雪夫多项式U0到U_{m_max}在x处的值
    U = zeros(m_max + 1, 1); % U(1)对应U0，U(2)对应U1，依此类推
    U(1) = 1; % U0(x) = 1
    if m_max >= 1
        U(2) = 2 * x; % U1(x) = 2x
        for m = 2:m_max
            % 递推公式: Um(x) = 2x*Um-1(x) - Um-2(x)
            U(m+1) = 2 * x * U(m) - U(m-1);
        end
    end
    
    % 计算各k >= n的n阶导数
    for k = n:(N-1)
        m = k - n;
        % 确保m在U数组的范围内
        if m > m_max
            continue; % 此情况理论上不会发生
        end
        
        % 计算下降阶乘(k)_n = k*(k-1)*...*(k-n+1)
        if k < n
            coeff = 0;
        else
            coeff = 1;
            for i = 0:(n-1)
                coeff = coeff * (k - i);
            end
        end
        
        % 计算导数: 2^{n-1} * (k)_n * U_{k-n}(x)
        dT(k+1) = (2^(n-1)) * coeff * U(m + 1);
    end
end
end